#include "modified_approach.h"
#include <iostream>
#include <device_launch_parameters.h>
#include <thrust/device_vector.h>
#include <thrust/sequence.h>
#include <thrust/sort.h>
#include <fstream>
#include <thrust/host_vector.h>
#include <thrust/scan.h>
#include <fstream>
#include <stdlib.h>     
#include <time.h> 
#include <math.h>
#include <set>
#include <sstream>
#include <device_launch_parameters.h>
#include <algorithm>
#include <cuda_profiler_api.h>

using namespace std;
using namespace thrust;


vector<double> exec_time;

#define MIN_VALUE -9999999.0;
int MAX_FEATURE_ID;

// keep track of the last positions the counting computed
struct QueryPos{
	int2 pos[MAX_DIM_X]; // pos[i].x means the down moved step, pos[i].y means the up moved step
	int2 detail_pos[MAX_DIM_X]; // detail position of 

	__host__ __device__ QueryPos(){
		for(int i = 0; i < MAX_DIM_X; i++){
			pos[i].x = 0;
			pos[i].y = 0;
			detail_pos[i].x = 0;
			detail_pos[i].y = 0;
		}
	}
};


// the entry of a inverted list
struct InvlistEnt{
	int feature_id; // this number must not exceed MAX_FEATURE_ID
	//float count;
};

// the entry of mapping from query to feature
struct QueryFeatureEnt{
	float count; // under current setup, count is actually integer and max is 128
	int ACD; // accumulated distance

	__host__ __device__ QueryFeatureEnt(){
		count = 0;
		ACD = 0;
	}
};

//Dec 26: read the whole inverted index from one binary file!
void get_invert_list_from_binary_file(string filename);
void init_problem_from_binary_file(string inv_list_file,string query_file_path, int num_of_query);

// random generate inverted list
void rand_inv_list();
// random generate query
void rand_query(int query_num);

// 
void get_invert_list_from_file(string data_root_dir);
void get_query_from_file(string query_file_name);

// print the problem
void print_problem();
void print_problem_host(
	host_vector<int>& query_keyword,
	host_vector<InvlistEnt>& invert_list,
	host_vector<int>& invert_list_idx);
void print_result(int qindex); // print query to feature mapping matrix
void print_final_result(vector<Result>& result);

bool compare_result(const Result& r1, const Result& r2){
	return r1.query < r2.query || (r1.query == r2.query && r1.feature_id < r2.feature_id);
}

// inverted list
device_vector<InvlistEnt>* d_invert_list; 
device_vector<int>* d_invert_list_idx;

// query to keyword: each query associate with 128 keywords strictly
device_vector<int>* d_query_keyword;

// result of query to feature result
device_vector<QueryFeatureEnt>* d_query_feature;

// positions of last counting
device_vector<QueryPos>* d_query_pos;

// record which query is still running, a entry is 1 if it is still running, 0 means the query has stopped
device_vector<int>* d_valid_query;
int last_valid_query; // all element before this index is still running NOT USED for FUTURE PURPOSE

int num_of_rounds;

__device__ __host__ void print_result_entry(Result& res){
	printf("%3d %5d %2d %9d %9d",res.query,res.feature_id,(int)res.count,res.lb,res.ub);
}

// distance function for values on the same dimension
template <class T>
__device__ T distance_func(T a, T b){
	int diff = a - b;
	return diff*diff;
}

// computing the counting while saving the position
__global__ void compute_mapping_saving_pos(
	int* query_keyword,
	InvlistEnt* invert_list, int* invert_list_idx,
	QueryFeatureEnt* query_feature,
	int N, // number of elements processd this round
	QueryPos* query_pos, // remeber the counting position
	int lb, // lower bound of the search, relative bounds 
	int ub, // upper bound of the search, relative bounds
	int MAX_FEATURE_ID
	){

	int bid = blockIdx.x; // correspond to query
	int tid = threadIdx.x;
	int block_size = blockDim.x;

	int start_idx = bid * MAX_DIM_X;

	// loading the position
	__shared__ int2 pos[MAX_DIM_X];
	__shared__ int2 detail_pos[MAX_DIM_X];
	if(tid < MAX_DIM_X){
		pos[tid] = query_pos[bid].pos[tid];
		detail_pos[tid] = query_pos[bid].detail_pos[tid];
	}

	__syncthreads();

	for(int i = 0; i < MAX_DIM_X; i++){
		int down_pos = pos[i].x;
		int up_pos = pos[i].y;
		int down_detail_pos = detail_pos[i].x;
		int up_detail_pos = detail_pos[i].y;
		
		int keyword = query_keyword[start_idx+i];
		int y_dim_value = keyword % MAX_DIM_Y;

		if(down_pos == 0 && up_pos == 0){
			int invert_list_start = keyword == 0 ? 0 : invert_list_idx[keyword-1];
			int invert_list_end = invert_list_idx[keyword];
			int invert_list_size = invert_list_end - invert_list_start;
			int process_round = invert_list_size / block_size + (invert_list_size % block_size != 0);

			for(int j = 0; j < process_round; j++){
				int idx = invert_list_start+j*block_size+tid;
				if(idx < invert_list_end){
					InvlistEnt inv_ent = invert_list[idx];
					int target_idx = bid*MAX_FEATURE_ID+inv_ent.feature_id;
					query_feature[target_idx].count += 1; //inv_ent.count;
					//printf("meet doc %d in point search in dim %d\n",inv_ent.feature_id,i);
					/*if(target_idx == 26*MAX_FEATURE_ID){
						printf("problem3 at tid %d and bid %d with count %d \n",tid,bid,inv_ent.count);
					}*/
				}
			}
			down_pos = 1;
			up_pos = 1;
		}


		// going downward
		int down_compute = N;
		int down_value = y_dim_value - down_pos; // move the position that last iteartion possessed.

		while(down_value >= 0 && down_compute > 0 && down_value > y_dim_value - lb){
			int down_keyword = down_value + i * MAX_DIM_Y;

			int invert_list_start = down_keyword == 0 ? 0 : invert_list_idx[down_keyword-1];
			invert_list_start += down_detail_pos;	// update starting position
			int invert_list_end = invert_list_idx[down_keyword];
			int remaining_size = invert_list_end - invert_list_start;
			
			int num_of_process_element = min(down_compute,remaining_size);
			int process_round = num_of_process_element / block_size + (num_of_process_element % block_size != 0);

			/*__syncthreads();
			if(threadIdx.x == 0) printf("\n num_of_process_element down is %d\n",num_of_process_element);
			__syncthreads();*/

			// calcuate the distance from the current point to the query along this dimension
			int dist_sqr = distance_func(y_dim_value,down_value);

			for(int j = 0; j < process_round; j++){
				int idx = invert_list_start+j*block_size+tid;
				if(idx < num_of_process_element + invert_list_start){
					InvlistEnt inv_ent = invert_list[idx];
					int target_idx = bid*MAX_FEATURE_ID+inv_ent.feature_id;
					query_feature[target_idx].count += 1; //inv_ent.count;
					query_feature[target_idx].ACD += dist_sqr;
					//printf("meet doc %d in down search in dim %d\n",inv_ent.feature_id,i);
					/*if(target_idx == 26*MAX_FEATURE_ID){
						printf("problem3 at tid %d and bid %d with count %d \n",tid,bid,inv_ent.count);
					}*/
				}
			}

			if(down_compute >= remaining_size){
				down_detail_pos = 0;
				down_compute -= num_of_process_element;
			}else{
				down_detail_pos += num_of_process_element;
				break;
			}
			down_value--;
		}

		// going upward
		int up_compute = N;
		int up_value = y_dim_value + up_pos; // move the position that last iteartion possessed.

		while(up_value < MAX_DIM_Y && up_compute > 0 && up_value < y_dim_value + ub){
			
			int up_keyword = up_value + i * MAX_DIM_Y;

			int invert_list_start = up_keyword == 0 ? 0 : invert_list_idx[up_keyword-1];
			invert_list_start += up_detail_pos;
			int invert_list_end = invert_list_idx[up_keyword];
			int remaining_size = invert_list_end - invert_list_start;
			
			int num_of_process_element = min(up_compute,remaining_size);
			int process_round = num_of_process_element / block_size + (num_of_process_element % block_size != 0);

			/*__syncthreads();
			if(threadIdx.x == 0) printf("\n num_of_process_element up is %d\n",num_of_process_element);
			__syncthreads();*/

			// calcuate the distance from the current point to the query along this dimension
			int dist_sqr = distance_func(up_value,y_dim_value);

			for(int j = 0; j < process_round; j++){
				int idx = invert_list_start+j*block_size+tid;
				if(idx < num_of_process_element + invert_list_start){
					InvlistEnt inv_ent = invert_list[idx];
					int target_idx = bid*MAX_FEATURE_ID+inv_ent.feature_id;
					query_feature[target_idx].count += 1; //inv_ent.count;
					query_feature[target_idx].ACD += dist_sqr;
					//printf("meet doc %d in up search in dim %d\n",inv_ent.feature_id,i);
					/*if(target_idx == 26*MAX_FEATURE_ID){
						printf("problem3 at tid %d and bid %d with count %d \n",tid,bid,inv_ent.count);
					}*/
				}
			}

			if(up_compute >= remaining_size){
				up_detail_pos = 0;
				up_compute -= num_of_process_element;
			}else{
				up_detail_pos += num_of_process_element;
				break;
			}

			up_value++;
		}

		//if(threadIdx.x == 0) printf("\n");
		__syncthreads();
		if(tid == 0){
			pos[i].x =  y_dim_value - down_value;
			pos[i].y = up_value - y_dim_value;
			detail_pos[i].x = down_detail_pos;
			detail_pos[i].y = up_detail_pos;
		}
	}

	if(tid < MAX_DIM_X){
		query_pos[bid].pos[tid] = pos[tid];
		query_pos[bid].detail_pos[tid] = detail_pos[tid];
	}

	// __syncthreads();
	//if(bid == 26 && tid == 0){
	//	printf("count 26 to 0 is %d\n",query_feature[26*MAX_FEATURE_ID].count);
	//	printf("MAX_FEATURE_ID %d\n",MAX_FEATURE_ID);
	//}
}

__global__ void prefix_count(
	QueryFeatureEnt* query_feature,
	int* count_vec, int* blk_count_vec, float threshold,int MAX_FEATURE_ID, int round_num){

	int bid = blockIdx.x;
	int tid = threadIdx.x;
	int oid = bid * blockDim.x + tid;

	int block_start = bid * MAX_FEATURE_ID;
	int block_end = block_start + MAX_FEATURE_ID;
	int round = MAX_FEATURE_ID / blockDim.x + (MAX_FEATURE_ID % blockDim.x != 0);


	__shared__ int shared_count;
	if(tid == 0) shared_count = 0;
	__syncthreads();

	int count = 0;
	for(int i = 0; i < round; i++){
		int idx = block_start + i * blockDim.x + tid;
		if(idx < block_end){
			if(query_feature[idx].count > threshold){
				count++;
				atomicAdd(&shared_count,1);
			}
		}
	}
	count_vec[oid] = count;
	__syncthreads();

	if(tid == 0){
		blk_count_vec[bid] = shared_count;
	}
}


__device__ void blk_sort_shared_memory(int* shared_array, int* temp_array, int size){
	int rounds = ceil(log2((double)size));

	int base = 1;
	for(int i = 1; i <= rounds; i++){
		
		if(threadIdx.x % (2*base) == 0){
			int start_idx_x = threadIdx.x;
			int end_idx_x = start_idx_x + base;
			end_idx_x = end_idx_x < size ? end_idx_x : size;

			int start_idx_y = end_idx_x;
			int end_idx_y = start_idx_y + base;
			end_idx_y = end_idx_y < size ? end_idx_y : size;

			int x_pointer = start_idx_x;
			int y_pointer = start_idx_y;

			int output_ptr = x_pointer;
			while(x_pointer < end_idx_x || y_pointer < end_idx_y){
				if(x_pointer >= end_idx_x)
					temp_array[output_ptr++] = shared_array[y_pointer++];
				else if(y_pointer >= end_idx_y)
					temp_array[output_ptr++] = shared_array[x_pointer++];
				else if(shared_array[x_pointer] < shared_array[y_pointer])
					temp_array[output_ptr++] = shared_array[y_pointer++];
				else
					temp_array[output_ptr++] = shared_array[x_pointer++];
			}			
		}
		__syncthreads();
		if(threadIdx.x < size){
			shared_array[threadIdx.x] = temp_array[threadIdx.x];
		}
		__syncthreads();
		base *= 2;
	}
}

// to test the block sorting of a shared memory
__global__ void test_shared_sort_device(int* result){
	__shared__ int input_array[100];
	__shared__ int temp_array[100];

	int start_idx = blockIdx.x * 100;
	if(threadIdx.x < 100) input_array[threadIdx.x] = result[start_idx + threadIdx.x];

	blk_sort_shared_memory(input_array,temp_array,100);

	if(threadIdx.x < 100) result[start_idx + threadIdx.x] = input_array[threadIdx.x];
}

// to test the block sorting of a shared memory
void test_shared_sort(){
	srand(2);
	device_vector<int> array_to_sort(1000);
	for(int i = 0; i < 1000; i++){
		array_to_sort [i] = rand() % 100;
	}

	for(int i = 0; i < 10; i++){
		for(int j = 0; j < 100; j++){
			cout << array_to_sort[i*100+j] << " ";
		}
		cout << endl;
	}

	cout << "----------------------------------" << endl;

	test_shared_sort_device<<<10,128>>>(raw_pointer_cast(array_to_sort.data()));
	cudaDeviceSynchronize();

	for(int i = 0; i < 10; i++){
		for(int j = 0; j < 100; j++){
			cout << array_to_sort[i*100+j] << " ";
		}
		cout << endl;
	}
}

__global__ void output_result(
	QueryFeatureEnt* query_feature,
	QueryPos* query_pos, // remeber the counting position
	int* query_keyword, // the information of query and it's value for each dimension
	int* ending_idx,float threshold,
	Result* result_vec,int MAX_FEATURE_ID, int round_num){

	int pos = blockIdx.x * blockDim.x + threadIdx.x;
	
	int start_idx = pos == 0 ? 0 : ending_idx[pos-1];

	int block_start = blockIdx.x * MAX_FEATURE_ID;
	int block_end = block_start + MAX_FEATURE_ID;
	int round = MAX_FEATURE_ID / blockDim.x + (MAX_FEATURE_ID % blockDim.x != 0);

	__shared__ int queryLowerBound[MAX_DIM_X];
	__shared__ int queryUpperBound[MAX_DIM_X];
	__shared__ int temp_shared[MAX_DIM_X];

	if(threadIdx.x < MAX_DIM_X){
		int2 q_pos = query_pos[blockIdx.x].pos[threadIdx.x];
		int max_lb = max(q_pos.x,q_pos.y);
		queryLowerBound[threadIdx.x] = max_lb * max_lb;

		int query_start_idx = blockIdx.x * MAX_DIM_X;
		int query_dim_value = (query_keyword[query_start_idx + threadIdx.x]) % MAX_DIM_Y;
		int max_ub =  max(query_dim_value,MAX_DIM_Y-query_dim_value);
		queryUpperBound[threadIdx.x] = max_ub * max_ub;
	}

	__syncthreads();

	// sort the lower and up bound
	blk_sort_shared_memory(queryLowerBound,temp_shared,MAX_DIM_X);
	blk_sort_shared_memory(queryUpperBound,temp_shared,MAX_DIM_X);

	//if(blockIdx.x == 26 && threadIdx.x == 0){
	//	for(int i = 0; i < MAX_DIM_X; i++){
	//		printf("%d ",queryLowerBound[i]);
	//	}
	//	printf("\n\n");
	//	for(int i = 0; i < MAX_DIM_X; i++){
	//		printf("%d ",queryUpperBound[i]);
	//	}
	//	printf("\n --------------------------");
	//}
	

	// calculate lower bound for this query, total distance for current retrieve point
	int lower_bound_sum = 0;
	for(int i = 0; i < MAX_DIM_X; i++) lower_bound_sum += queryLowerBound[i];

	for(int i = 0; i < round; i++){
		int idx = block_start + i * blockDim.x + threadIdx.x;
		if(idx < block_end){
			float count = query_feature[idx].count;
			if(count > threshold){
				Result new_result;
				new_result.query = blockIdx.x;
				new_result.feature_id = i * blockDim.x + threadIdx.x;
				new_result.count = count;
				
				int ACD = query_feature[idx].ACD;

				// DL(f) = DL(us) - sum(max count number of dt) + ACD(f)
				int lb = lower_bound_sum;
				for(int m = 0; m < (int)count; m++){
					lb -= queryLowerBound[m];
				}
				lb += ACD;

				// DU(f) = ACD(f) + sum(max 128 - c of upper bound)
				int ub = ACD;
				for(int m = 0; m < MAX_DIM_X - (int)count; m++){
					ub += queryUpperBound[m];
				}

				new_result.lb = lb;
				new_result.ub = ub;

				result_vec[start_idx] = new_result;
				start_idx++;
			}
		}
	}
}

__global__ void terminate_check(Result* ub_sorted,Result* lb_sorted, int* end_idx, int* output, int K, int round_num){
	int tid = threadIdx.x;
	int bid = blockIdx.x;

	int blk_start_idx = bid == 0 ? 0 : end_idx[bid-1];
	int blk_end_idx = end_idx[bid];

	/*if(tid == 0 && bid == 26){
		printf("reached here: bid: %d with block start idx %d and end idx %d \n",bid,blk_start_idx,blk_end_idx);
	}*/

	if(blk_end_idx - blk_start_idx <= K){
		return;
	}
/*
	if(round_num == 18 && bid == 26 && tid == 0){
		for(int i = 0; i < K+1; i++){
			print_result_entry(ub_sorted[blk_start_idx+i]);
			printf("\n");
		}
		printf("###############################\n");
		for(int i = 0; i < K+1; i++){
			print_result_entry(lb_sorted[blk_start_idx+i]);
			printf("\n");
		}
	}
*/
	blk_end_idx = blk_start_idx + K + 1;

	

	int round = (K+1) / blockDim.x + ((K+1) % blockDim.x != 0);

	__shared__ int termination_check;
	if(tid == 0) termination_check = 0;
	__syncthreads();

	for(int i = 0; i < round; i++){
		int idx = blk_start_idx + i * blockDim.x + tid;
		if(idx < blk_end_idx){
			Result lb_elem = lb_sorted[idx];
			bool duplicated = false; // if element in the lb sorted array existed in the ub sorted array
			for(int j = 0; j < K; j++){
				Result ub_elem = ub_sorted[blk_start_idx + j];
				if(lb_elem.feature_id == ub_elem.feature_id){
					duplicated = true;
					break;
				}
			}
			Result ub_elem = ub_sorted[blk_start_idx + K - 1];
			bool lb_greater_ub = (lb_elem.lb >= ub_elem.ub);

			if(round_num == 18 && bid == 26){
				printf("tid %d has ub_elem: %d and lb_elem: %d with ub: %d and lb: %d and duplicated: %d\n",
					tid,ub_elem.feature_id,lb_elem.feature_id,ub_elem.ub,lb_elem.lb,duplicated);
			}

			if(!duplicated && lb_greater_ub){
				//printf("reached here tid %d and bid %d\n",tid,bid);
				atomicAdd(&termination_check,1);
	
			
			}
		}
	}
	__syncthreads();
	if(tid == 0){
		if(termination_check != 0){
			output[bid] = 0;
			
		}
	}
}

__global__ void extract_topK(Result* ub_sorted,int* end_idx,Result* output, int K){
	int tid = threadIdx.x;
	int bid = blockIdx.x;
	int oid = bid * K + tid;

	int blk_start_idx = bid == 0 ? 0 : end_idx[bid-1];
	int blk_end_idx = blk_start_idx + K;

	int round = K / blockDim.x + (K % blockDim.x != 0);

	for(int i = 0; i < round; i++){
		int idx = blk_start_idx + i * blockDim.x + tid;
		if(idx < blk_end_idx){
			output[oid] = ub_sorted[idx];
			oid += blockDim.x;
		}
	}
}

double get_uniform_rand(){
	return (double)(rand() % 100000) / 100000;
}

template <class T>
void init_gpu_memory(host_vector<T>& hvec, device_vector<T>** dvec){
	*dvec = new device_vector<T>(hvec.size());
	**dvec = hvec;
}

void rand_inv_list(){	
	MAX_FEATURE_ID = 100;
	vector<vector<int> > keyword_feature(MAX_DIM_X*MAX_DIM_Y); 
	for(int i = 0; i < MAX_FEATURE_ID; i++){
		for(int x = 0; x < MAX_DIM_X;  x++){
			int keyword = x * MAX_DIM_Y + rand() % MAX_DIM_Y;
			keyword_feature[keyword].push_back(i);
		}
	}

	host_vector<InvlistEnt> h_invert_list;
	host_vector<int> h_invert_list_idx;
	int count = 0;
	// construct inverted list
	for(int i = 0; i < MAX_DIM_X*MAX_DIM_Y; i++){
		for(int j = 0; j < keyword_feature[i].size(); j++){
			InvlistEnt new_ent;
			new_ent.feature_id = keyword_feature[i][j];
			//new_ent.count = 1;
			h_invert_list.push_back(new_ent);
		}
		count += keyword_feature[i].size();
		h_invert_list_idx.push_back(count);
	}
	
	// transfer to gpu
	init_gpu_memory(h_invert_list,&d_invert_list);
	init_gpu_memory(h_invert_list_idx,&d_invert_list_idx);
}

void get_query_from_file(string query_file_name){
	ifstream query_file(query_file_name.c_str());
	if(!query_file.is_open()){
		cout << "query file cannot found" << endl;
		exit(1);
	}

	host_vector<int> h_query_keyword;
	int num_query = 0;
	while(!query_file.eof()){
		string input_line;
		getline(query_file,input_line);
		stringstream ss_test(input_line);

		string test;
		ss_test >> test;
		if(test.empty()) continue;

		num_query++;
		stringstream ss(input_line);

		for(int i = 0; i < MAX_DIM_X; i++){
			int keyword = -1;
			ss >> keyword;
			h_query_keyword.push_back(keyword);
		}
	}

	d_query_pos = new device_vector<QueryPos>(num_query);
	init_gpu_memory(h_query_keyword,&d_query_keyword);
	cout << "finish reading query from file" << endl;
}	

void rand_query(int query_num){

	host_vector<int> h_query_keyword;
	// construct feature to keyword
	for(int i = 0; i < query_num; i++){
		for(int x = 0; x < MAX_DIM_X; x++){
			int new_keyword = x * MAX_DIM_Y + rand() % MAX_DIM_Y;
			h_query_keyword.push_back(new_keyword);
		}
	}
	
	d_query_pos = new device_vector<QueryPos>(query_num);

	init_gpu_memory(h_query_keyword,&d_query_keyword);
	//cout << "finish rand query" << endl;
}


//Dec 27: read inverted lists from a single binary file
void init_problem_from_binary_file(string inv_list_file,string query_file_path, int num_of_query){

	if(inv_list_file.empty()) rand_inv_list();
	else get_invert_list_from_binary_file(inv_list_file);

	if(query_file_path.empty()) rand_query(num_of_query);
	else get_query_from_file(query_file_path);

	
	
	double total_bytes = (double)d_invert_list->size() * sizeof(InvlistEnt) + 
		d_invert_list_idx->size() * sizeof(int) + 
		d_query_keyword->size() * sizeof(int) + 
		d_query_pos->size() * sizeof(QueryPos) + 
		2*MAX_FEATURE_ID*d_query_pos->size()*sizeof(int);

        cout << "total_memory used: " << total_bytes / pow(2.0,30.0) << "GB" << endl;

      cout << "inverted index: " << (double)d_invert_list->size() * sizeof(InvlistEnt) /pow(2.0,30.0) << " GB" <<endl;

	if((double)MAX_FEATURE_ID*d_query_pos->size()*sizeof(int) > 5.0 * pow(2.0,30.0)){
		cout << "too much memory is used for query feature mapping" << endl;
		exit(1);
	}
	//d_query_feature = new device_vector<float>((unsigned int)MAX_FEATURE_ID*d_query_pos->size(),0);
	d_query_feature = new device_vector<QueryFeatureEnt>((unsigned int)MAX_FEATURE_ID*d_query_pos->size());	
	d_valid_query = new device_vector<int>(d_query_pos->size(),1);
	num_of_rounds = 0;
	
}

//Dec 27
void get_invert_list_from_binary_file(string filename)
{
	
	// initialize the previousKey to be 0;
		int previousKey = 0;
		int key;
		int numberOfFeatures;
		int keyPosition = 0;		
		MAX_FEATURE_ID = -1;			
				
		vector<InvlistEnt>  featureIdList;
		vector<int>  featureIdListIndex;				
		
		clock_t start = clock();							  
		ifstream inFile(filename.c_str(), ios::in | ios::binary );
		if ( ! inFile )
		{
			cerr << "error : unable to open input file: " << inFile << endl;
		}
		else if ( ! inFile.is_open() )
		{
			cerr << "error : unable to open input file: " << inFile << endl;
		}
		else
		{
			// read the key, if reach the end of the file, jump out
			while ( inFile.read((char*) (&key), sizeof(int) ) )
			{				
				// add the position for empty keys and the current key
				for ( int j = previousKey; j < key; j++ )
				{
					featureIdListIndex.push_back(keyPosition);
				}

				// read the number of features
				inFile.read((char*) (&numberOfFeatures), sizeof(int) );
				
				// update the current position, and the previousKey
				keyPosition += numberOfFeatures;
				previousKey = key;

				// read the feature ids, and fill in the featureIdList
				if ( numberOfFeatures > 0 )
				{
					int featureID;			
					for ( int j = 0; j < numberOfFeatures; j++ )
					{
						// read the feature ID list
						inFile.read((char*) (&featureID), sizeof(int) );
						
						//if(featureID == 0)
						//	cout << "key: " << key << ", ";
						
						InvlistEnt new_ent;
						new_ent.feature_id = featureID;
						if(new_ent.feature_id > MAX_FEATURE_ID){
							MAX_FEATURE_ID = new_ent.feature_id;
						}
						featureIdList.push_back(new_ent);
					}
				}								
			}
			
			MAX_FEATURE_ID++;  //very important!
			
			// fill the key positions for the empty keys
			
			for ( int j = previousKey ; j < MAX_DIM_X * MAX_DIM_Y; j++)
			{
				featureIdListIndex.push_back(keyPosition);
			}

			inFile.close();
			inFile.clear();
			
		}
		
		double init_time = (double)(clock() - start) / CLOCKS_PER_SEC;
		cout<< "Reading Inverted Lists from disk to Host Memory takes: " << init_time << " seconds. " <<endl;
				
		start = clock();
		host_vector<InvlistEnt> featureIdList_h(featureIdList); 
		host_vector<int>  featureIdListIndex_h(featureIdListIndex);
		cout << "size of list: " << featureIdList_h.size() * sizeof(InvlistEnt) << endl;
		cout<< "host vector allocated" <<endl;
		init_gpu_memory(featureIdList_h,&d_invert_list);
		init_gpu_memory(featureIdListIndex_h,&d_invert_list_idx);
		init_time = (double)(clock() - start) / CLOCKS_PER_SEC;
		cout<< "Loading from CPU to GPU takes: " << init_time << " seconds. " <<endl;
	
}

void get_invert_list_from_file(string data_root_dir){

	string sorted_file_name(data_root_dir);
	sorted_file_name.append("sorted.txt");
	ifstream sorted_file(sorted_file_name.c_str());
	if(!sorted_file.is_open()){
		cout << "sorted index file not found" << endl;
		exit(1);
	}

	host_vector<InvlistEnt> h_inv_list;
	host_vector<int> h_inv_list_idx;
	int count = 0;
	int keyword_idx = 0;

	MAX_FEATURE_ID = -1;

	while(!sorted_file.eof()){
		string input_line;
		getline(sorted_file,input_line);
		stringstream ss(input_line);

		string name;
		int file_idx;
		int file_len;

		ss >> name;
		if(name.empty()) continue;
		ss >> file_len;

		file_idx = atoi(name.c_str());
		string file_name(data_root_dir);
		file_name.append("sift_small/");
		file_name.append(name.c_str());

		for(int i = keyword_idx; i < file_idx; i++){
			h_inv_list_idx.push_back(count);
		}

		keyword_idx = file_idx;
		count += file_len;

		ifstream keyword_file(file_name.c_str());
		if(!keyword_file.is_open()){
			cout << "invert list file not found" << endl;
			exit(1);
		}

		for(int i = 0; i < file_len; i++){
			InvlistEnt new_ent;
			keyword_file >> new_ent.feature_id;
			//new_ent.count = 1;
			h_inv_list.push_back(new_ent);
			if(new_ent.feature_id > MAX_FEATURE_ID){
				MAX_FEATURE_ID = new_ent.feature_id;
			}
		}

		keyword_file.close();

	} 
	//cout << "finish reading file" << endl;
	//cout << "max feature id: " << MAX_FEATURE_ID << endl;
	sorted_file.close();

	MAX_FEATURE_ID++;

	for(int i = keyword_idx; i < MAX_DIM_X * MAX_DIM_Y; i++){
		h_inv_list_idx.push_back(count);
	}

	init_gpu_memory(h_inv_list,&d_invert_list);
	init_gpu_memory(h_inv_list_idx,&d_invert_list_idx);
}

void init_problem(string inv_list_dir,string query_file_path, int num_of_query){

	if(inv_list_dir.empty()) rand_inv_list();
	else get_invert_list_from_file(inv_list_dir);

	if(query_file_path.empty()) rand_query(num_of_query);
	else get_query_from_file(query_file_path);

	/*double total_bytes = (double)d_invert_list->size() * sizeof(InvlistEnt) + 
		d_invert_list_idx->size() * sizeof(int) + 
		d_query_keyword->size() * sizeof(int) + 
		d_query_pos->size() * sizeof(QueryPos) + 
		MAX_FEATURE_ID*d_query_pos->size()*sizeof(int);*/

	if((double)MAX_FEATURE_ID*d_query_pos->size()*sizeof(QueryFeatureEnt) > 5.0 * pow(2.0,30.0)){
		cout << "too much memory is used for query feature mapping" << endl;
		exit(1);
	}
	d_query_feature = new device_vector<QueryFeatureEnt>((unsigned int)MAX_FEATURE_ID*d_query_pos->size());
	d_valid_query = new device_vector<int>(d_query_pos->size(),1);

	num_of_rounds = 0;
}

bool compute_feature_mapping(vector<Result>& result,int N,int lb,int ub,float threshold, int K){
	if(THREAD_PER_BLK < MAX_DIM_X){
		cout << "thread per block must be larger than max dimension of x" << endl;
		exit(1);
	}
	cudaProfilerStart();
	cout << "started counting" << endl;
	cout << "max feature ID: " << MAX_FEATURE_ID << endl;

	cudaEvent_t start, stop;
	float elapsedTime;
	cudaEventCreate(&start);
	cudaEventCreate(&stop);

	cudaEventRecord(start, 0);

	compute_mapping_saving_pos<<<d_query_pos->size(),THREAD_PER_BLK>>>(
		raw_pointer_cast(d_query_keyword->data()),
		raw_pointer_cast(d_invert_list->data()),
		raw_pointer_cast(d_invert_list_idx->data()),
		raw_pointer_cast(d_query_feature->data()),
		N,raw_pointer_cast(d_query_pos->data()),lb,ub,MAX_FEATURE_ID);

		//thrust::host_vector<QueryFeatureEnt> dvec(d_query_feature->begin(),d_query_feature->end());
		/*thrust::host_vector<InvlistEnt> dinv(d_invert_list->begin(),d_invert_list->end());
		int count_f = 0;
		for(int i = 0; i < dinv.size(); i++){
			if(dinv[i].feature_id == 0){
				count_f += dinv[i].count;
			}
		}
		cout << "count_f is " << count_f << endl;
		QueryFeatureEnt ent = (*d_query_feature)[26*MAX_FEATURE_ID];
		cout << "count at query 26 and feature 0:" << ent.count << endl;
		cout << "%%%%%%%%%%%%%%%%" << endl;*/

	device_vector<int> threshold_count(d_query_pos->size()*THREAD_PER_BLK,0);
	device_vector<int> query_result_count(d_query_pos->size(),0);
	cout << "started prefix counting" << endl;

	prefix_count<<<d_query_pos->size(),THREAD_PER_BLK>>>(
		raw_pointer_cast(d_query_feature->data()),
		raw_pointer_cast(threshold_count.data()),
		raw_pointer_cast(query_result_count.data()),
		threshold,MAX_FEATURE_ID,num_of_rounds);

	//if(num_of_rounds == 18){
		/*thrust::host_vector<QueryFeatureEnt> dvec(d_query_feature->begin(),d_query_feature->end());
		cout << "count at query 26 and feature 0:" << dvec[26*MAX_FEATURE_ID].count << endl;
		cout << "%%%%%%%%%%%%%%%%" << endl;*/
	//}
		/*thrust::host_vector<QueryPos> ddddvec(d_query_pos->begin(),d_query_pos->end());
		for(int i = 0; i < 128; i++){
			cout << "(" << ddddvec[26].pos[i].x << "," << ddddvec[26].pos[i].y << ") ";
		}
		cout << endl;
		*/
	//	thrust::host_vector<InvlistEnt> dinv(d_invert_list->begin(),d_invert_list->end());
		/*thrust::host_vector<int> dinvidx(d_invert_list_idx->begin(),d_invert_list_idx->end());
		cout << "invert list size: " << dinvidx.size() << endl;*/

		/*for(int i = 0; i < 128; i++){
			int count = 0;
			for(int j = 0; j < 256; j++){
				int key = i*256+j;
				int start_idx = (key == 0) ? 0 : dinvidx[key-1];
				int end_idx = dinvidx[key];
				for(int k = start_idx; k < end_idx; k++){
					if(dinv[k].feature_id == 0) count++;
				}
			}
			if(count != 1) cout << " dim: " << i << " has zero: " << count;
		}*/
		
	thrust::inclusive_scan(threshold_count.begin(),threshold_count.end(),threshold_count.begin()); // per thread inclusive scan
	thrust::inclusive_scan(query_result_count.begin(),query_result_count.end(),query_result_count.begin()); // per block inclusive scan

	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop); 
	exec_time[0] += elapsedTime;
	
	
	cout << "allocate d_result_ub_sorted" << endl;
	device_vector<Result> d_result_ub_sorted(threshold_count[threshold_count.size()-1]);
	
	cout << "start output result" << endl;

	
	cudaEventRecord(start, 0);
	
	output_result<<<d_query_pos->size(),THREAD_PER_BLK>>>(
		raw_pointer_cast(d_query_feature->data()),
		raw_pointer_cast(d_query_pos->data()),
		raw_pointer_cast(d_query_keyword->data()),
		raw_pointer_cast(threshold_count.data()),
		threshold,
		raw_pointer_cast(d_result_ub_sorted.data()),MAX_FEATURE_ID,num_of_rounds);
	
	
	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop); 
	exec_time[1] += elapsedTime;
	
	cout << "allocate d_result_lb_sorted" << endl;
	device_vector<Result> d_result_lb_sorted(d_result_ub_sorted.begin(),d_result_ub_sorted.end());

	/*cout << "number of total result: " << d_result_ub_sorted.size() << endl;
	clock_t start = clock();*/

	cudaEventRecord(start, 0);

	thrust::sort(d_result_ub_sorted.begin(),d_result_ub_sorted.end(),Ubcomapre());
	thrust::sort(d_result_lb_sorted.begin(),d_result_lb_sorted.end(),Lbcompare());

	thrust::host_vector<Result> h_result = d_result_ub_sorted;
	for(int i = 0; i < 5; i++){
		 print_result_entry(h_result[i]);
		 cout << endl;
	}
	//host_vector<int> keyword_vec = *d_query_keyword;
	//host_vector<QueryPos> pos_vec = *d_query_pos;
	//for(int i = 0 ; i < MAX_DIM_X ; i++){
	//	cout << "dim " << i << "has value " << (keyword_vec[i]%MAX_DIM_Y) << " : " << pos_vec[0].pos[i].x << "(" << pos_vec[0].detail_pos[i].x << ") --- " << pos_vec[0].pos[i].y << "(" << pos_vec[0].detail_pos[i].y << ")" << endl;
	//}
	cout << "+++++++++++++++++++++++" << endl;

	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop); 
	exec_time[2] += elapsedTime;
	cout << "finish sort" << endl;
	/*double time = (double)(clock() - start) / CLOCKS_PER_SEC;
	cout << "sort timing: " << time << endl;*/

	//for(int i = 0; i < d_result_ub_sorted.size(); i++){
	//	Result ub_res = d_result_ub_sorted[i];
	//	Result lb_res = d_result_lb_sorted[i];
	//	print_result_entry(ub_res);
	//	cout << " | ";
	//	print_result_entry(lb_res);
	//	cout << endl;
	//}
	cudaEventRecord(start, 0);
	terminate_check<<<d_query_pos->size(),THREAD_PER_BLK>>>(
		raw_pointer_cast(d_result_ub_sorted.data()),
		raw_pointer_cast(d_result_lb_sorted.data()),
		raw_pointer_cast(query_result_count.data()),
		raw_pointer_cast(d_valid_query->data()),K,num_of_rounds);  

	int terminate_sum = thrust::reduce(d_valid_query->begin(),d_valid_query->end());
	cout << "finish reduce" << endl;
	
	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop); 
	exec_time[3] += elapsedTime;
	cudaEventDestroy(start);
	cudaEventDestroy(stop);

	cout << "number of total result: " << d_result_ub_sorted.size() << " with " << terminate_sum << " queries still running in rounds " << num_of_rounds << endl;

	cudaProfilerStop();
	num_of_rounds++;
	

	if(terminate_sum == 0){ // means all queries are terminated
		device_vector<Result> d_result(d_query_pos->size()*K);
		extract_topK<<<d_query_pos->size(),THREAD_PER_BLK>>>(
			raw_pointer_cast(d_result_ub_sorted.data()),
			raw_pointer_cast(query_result_count.data()),
			raw_pointer_cast(d_result.data()),K);

		//return result
		result.resize(d_result.size());
		cudaMemcpy(&result[0],raw_pointer_cast(d_result.data()),d_result.size()*sizeof(Result),cudaMemcpyDeviceToHost);
		return true;
	}

	return false;
	//result.resize(d_result.size());
	//cudaMemcpy(&result[0],raw_pointer_cast(d_result.data()),d_result.size()*sizeof(Result),cudaMemcpyDeviceToHost);

}

void clear_gpu(){
	delete d_invert_list;
	delete d_invert_list_idx;
	delete d_query_keyword;
	delete d_query_feature;
	delete d_query_pos;
	delete d_valid_query;
}

void print_problem(){
	cout << "query to keyword relationship" << endl;
	for(int i = 0; i < d_query_pos->size(); i++){
		cout << "Q " << i << ":";
		for(int j = 0; j < MAX_DIM_X; j++){
			int keyword = (*d_query_keyword)[i*MAX_DIM_X+j];
			int x = keyword / MAX_DIM_Y;
			int y = keyword % MAX_DIM_Y;
			cout << "(" << x << "," << y << ");";
		}
		cout << endl;
	}
	cout << "inverted list" << endl;
	for(int x = 0; x < MAX_DIM_X; x++){
		for(int y = 0; y < MAX_DIM_Y; y++){
			

			int keyword = x*MAX_DIM_Y+y;
			int start_idx = (keyword == 0) ? 0 : (*d_invert_list_idx)[keyword-1];
			int end_idx = (*d_invert_list_idx)[keyword];

			if(start_idx == end_idx) continue;

			cout << "(" << x << "," << y << "):";
			for(int i = start_idx; i < end_idx; i++){
				InvlistEnt ent = (*d_invert_list)[i];
				cout << ent.feature_id << " ";
			}
			cout << endl;
		}
	}
}

void print_problem_host(
	host_vector<int>& query_keyword,
	host_vector<InvlistEnt>& invert_list,
	host_vector<int>& invert_list_idx){

	cout << "query to keyword relationship" << endl;
	for(int i = 0; i < d_query_pos->size(); i++){
		cout << "Q " << i << ":";
		for(int j = 0; j < MAX_DIM_X; j++){
			int keyword = query_keyword[i*MAX_DIM_X+j];
			int x = keyword / MAX_DIM_Y;
			int y = keyword % MAX_DIM_Y;
			cout << "(" << x << "," << y << ");";
		}
		cout << endl;
	}
	cout << "inverted list" << endl;
	for(int x = 0; x < MAX_DIM_X; x++){
		for(int y = 0; y < MAX_DIM_Y; y++){
	
			int keyword = x*MAX_DIM_Y+y;
			int start_idx = (keyword == 0) ? 0 : invert_list_idx[keyword-1];
			int end_idx = invert_list_idx[keyword];

			if(start_idx == end_idx) continue;

			cout << "(" << x << "," << y << "):";

			for(int i = start_idx; i < end_idx; i++){
				InvlistEnt ent = invert_list[i];
				cout << ent.feature_id << " ";
			}
			cout << endl;
		}
	}
}

void print_result(int qindex){
	if(qindex < 0){
		cout << "result for all" << endl;
		for(int i = 0; i < d_query_pos->size(); i++){
			for(int j =0; j < MAX_FEATURE_ID; j++){
				QueryFeatureEnt result = (*d_query_feature)[i*MAX_FEATURE_ID+j];
				cout << "(" << result.count << "," << result.ACD << " ";
			}
			cout << endl;
		}
	}else{
		cout << "result for query " << qindex << endl;
		for(int j =0; j < MAX_FEATURE_ID; j++){
			QueryFeatureEnt result = (*d_query_feature)[qindex*MAX_FEATURE_ID+j];
			cout << "(" << result.count << "," << result.ACD << " ";
		}
		cout << endl;
	}
}

void print_final_result(vector<Result>& result){
	cout << "result size is: " << result.size() << endl;
	std::sort(result.begin(),result.end(),compare_result);
	int qindex = -1;
	for(int i = 0; i < result.size(); i++){
		Result new_result = result[i];
		if(new_result.query != qindex){
			cout << endl;
			cout << new_result.query << ":" << "(" << new_result.feature_id << "," << new_result.count << ") ";
			qindex = new_result.query;
		}else{
			cout << "(" << new_result.feature_id << "," << new_result.count << ") ";
		}
	}
}

int main(){
	srand(1);

	vector<Result> result;
	int N = 10000; //1000;
	int num_of_query = 10; //100;
	int lb = 256;
	int ub = 256;
	float threshold = 0;
	int K = 3;
	
	clock_t start = clock();
	//init_problem("","",1);
	//init_problem("../data/sift_exp_data/","",num_of_query);
	init_problem_from_binary_file("../data/sift_4m/inverted_index.bin","query.txt",num_of_query);
	//init_problem_from_binary_file("../data/sift_1m/1m_inverted.bin","",num_of_query);
	double init_time = (double)(clock() - start) / CLOCKS_PER_SEC;
    cout<< "Initilization takes: " << init_time << " seconds. " <<endl;
	//cout << "finish init" << endl;


	//host_vector<InvlistEnt> h_inv_list = *d_invert_list;
	//host_vector<int> h_inv_list_idx = *d_invert_list_idx;

	//for(int i = 0; i < MAX_DIM_X; i++){
	//	for(int j = 0; j < MAX_DIM_Y; j++){
	//		int keyword = i*MAX_DIM_Y+j;
	//		int start = (keyword == 0) ? 0 : h_inv_list_idx[keyword-1];
	//		int end = h_inv_list_idx[keyword];
	//		cout << "(" << i << "," << j << "):";
	//		for(int k = start; k < end; k++) cout << h_inv_list[k].feature_id << " ";
	//		cout << endl;
	//	}
	//	cout << endl;
	//}

	exec_time.resize(4,0);
	
	double total_time = 0;
	bool finished = false;
	while(!finished){
		
		clock_t start = clock();
		finished = compute_feature_mapping(result,N,lb,ub,threshold,K);
		cudaDeviceSynchronize();
		double time = (double)(clock() - start) / CLOCKS_PER_SEC;
		total_time += time;
		//getchar();

		cout << "iteration takes time: " << time << endl;
	}

	cout << "finished with total time : " << total_time << endl;

	for(int i = 0; i < result.size(); i++){
		Result res = result[i];
		print_result_entry(res);
		cout << endl;
	}

	cout << "timing" << endl;
	for(int i = 0; i < exec_time.size(); i++){
		cout << "part " << i << " takes " << exec_time[i] << endl;
	}

	//init_problem("../data/sift_exp_data/","../data/sift_exp_data/query_format1/002_0001.sift",0);
	//cudaProfilerStart();

	//cudaProfilerStop();
	//print_final_result(result1);
	
	//print_final_result(result2);

	
	return 1;
}
